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We present a general framework to describe the evolutionary dynamics of an arbitrary number of types in 
finite populations based on stochastic differential equations (SDE). For large, but finite populations this allows 
to include demographic noise without requiring explicit simulations. Instead, the population size only rescales 
the amplitude of the noise. Moreover, this framework admits the inclusion of mutations between different types, 
provided that mutation rates, /i, are not too small compared to the inverse population size 1/A'^. This ensures 
that all types are almost always represented in the population and that the occasional extinction of one type does 
not result in an extended absence of that type. For fiN <C 1 this limits the use of SDE's, but in this case there 
are well established alternative approximations based on time scale separation. We illustrate our approach by 
a Rock-Scissors-Paper game with mutations, where we demonstrate excellent agreement with simulation based 
results for sufficiently large populations. In the absence of mutations the excellent agreement extends to small 
population sizes. 
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PACS numbers: 87.23.-n, 89.65.-s 02.50.Ey 
I. INTRODUCTION 

Populations evolve when different individuals have differ- 
ent traits or strategies that determine their reproductive suc- 
cess. In population genetic models, this success is typically 
constant, whereas in evolutionary game theory the fitness of 
an individual depends on interactions with other members of 
the population. Consequently, fitness depends on the rela- 
tive proportions (or frequencies) of different strategic types 
and hence gives rise to so called frequency dependent selec- 
tion. More successful strategies increase in abundance and 
may either take over the entire population or, as the strategy 
becomes increasingly common, may suffer from a decrease in 
fitness, which can result in the co-existence of two or more 
traits within the population. For example, in host-parasite co- 
evolution, a rare parasite may be most successful, but once 
it reaches high abundance, selection pressure on the host in- 
creases and the host is likely to develop some defense mech- 
anism. Consequently, the success of the parasite decreases. 
Such dynamics can be conveniently described by evolution- 
ary game theory i^^. 

Traditionally, the mathematical description of evolutionary 
game dynamics is formulated in terms of the deterministic 
replicator equation [2]. This implies that population sizes 
are infinite and populations are unstructured. Only more 
recently the stochastic dynamics in finite populations at- 
tracted increasing attention ITllil]. Typically, the stochastic 
approach becomes deterministic in the limit of infinite 
populations. For large, but finite populations the dynamics 
can be approximated by stochastic differential equations 
(SDE) Isl- UOll . This analytic approach is a natural extension 
of the replicator dynamics, which is capable of bridging 
the gap between deterministic models and individual based 
simulations. Moreover, SDE's are typically computationally 



far less expensive than simulations because the execution 
time does not scale with population size. 

II. MASTER EQUATION 

In unstructured, finite populations of constant size, N, con- 
sisting of d distinct strategic types and with a mutation rate, /i, 
evolutionary changes can be described by the following class 
of birth-death processes: In each time step, one individual of 
type j produces a single offspring and displaces another ran- 
domly selected individual of type k. With probability 1 — /i, 
no mutation occurs and j produces an offspring of the same 
type. But with probability /i, the offspring of an individual 
of type i {i ^ j) mutates into a type j individual. This re- 
sults in two distinct ways to increase the number of j types 
by one at the expense of decreasing the number of k types by 
one, hence keeping the population size constant. Biologically, 
keeping N constant implies that the population has reached 
a stable ecological equilibrium and assumes that this equilib- 
rium remains unaffected by trait frequencies. The probability 
for the event of replacing a type k individual with a type j 
individual is denoted by Tkj and is a function of the state of 
the population X = (Xi, X2, . . . Xd), with X„ indicating the 
number of individuals of type n such that X]n=i -^n ~ N. 

For this process it is straight forward to write down a 
Master equation iflOll and, at this point, there is no need to 
further specify the transition probabilities T^j. 

d 

pr+\x)=p\x)+Y,P^{x';)Tkj{x';)-p-{x)T,k{x)) 
j,k=i 

where P'^{X) denotes the probability of being in state JO at 
time r and X^ == {Xi , . . . X^ - 1, . . . X^ + 1, . . . Xd) repre- 
sents a state adjacent to X. 



III. FOKKER-PLANCK AND LANGEVIN EQUATIONS 

While the Master equation ([T]i is rather unwieldy, a 
Kramers-Moyal expansion yields a convenient approximation 
for large but finite N in the form of a Fokker-Planck equation 
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-p{x)Bjkix) 
(2) 



where x ~ X /N represents the state of the population in 
terms of frequencies of the different strategic types and p{x) 
is the probability density in state x. Due to the normalization 
Sfe=i ^k = 1, it suffices to consider d — 1 elements of the 
deterministic drift vector Akix), fc = 1, . . . d — 1. Similarly, 
we only need to consider a diffusion matrix Bjk{x) with di- 
mension (d — 1) X (d — 1), i.e. j,k = 1, . . .d — 1. The drift 
vector Ak{x) is given by 



Ak{x) 
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E 



T,k{x)-Tkj{x) 






(3) 



For the second equality we have used J2j=i Tkj{x) = 1, 
which simply states that a fc-type individual transitions to 
some other type (including staying type k) with probability 
one. Ak{x) is bounded in [—l,d — 1] because the Tjk are 
probabilities. 

The diffusion matrix Bjk{x) is defined as 

Bjk{x) = -^[T,k{x)+Tkjix)] for j^k (4a) 



Bjj{x) 
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J2 [T,i{x)+Ti,{x) 
1=1.1^ J 



(4b) 



For a detailed derivation see e.g. lIToll or 1IT2I1 . Note that the 
diffusion matrix is symmetric, Bjk{x) = Bkj{x) and van- 
ishes as ~ 1/A^ in the limit N ^ 00. Moreover, we have for 
all diagonal elements Bjj (x) > and for the non-diagonal el- 
ements Bjk{x) < 0. Note that Tjjix) cancels in Eq. (|3]l and 
does not appear in Eq. (|4|i - it is thus of no further concern. 

The noise of our underlyingprocess is uncorrected in time 
and hence the Ito calculus llllfl can be applied to derive a 
Langevin equation, which represents in our case a stochastic 
replicator-mutator equation. 



transition probabilities the matrix C{x) provides a quantita- 
tive description of the fluctuations introduced by microscopic 
processes in finite populations. In the limit N —>■ 00 the ma- 
trix C{x) vanishes with ~ I/VN and we recover a deter- 
ministic replicator mutator equation. Note that the replicator 
equation does not impose an upper or lower bound on A (c.f. 
Eq. ([3])). However, this difference merely amounts to a (con- 
stant) rescaling of time. 

The multiplicative character of the noise and its strategy- 
strategy correlations are determined by the form of the ma- 
trix C{x). In order to determine C{x), we first diagonalize 
B{x). Because B{x) is real and symmetric it is diagonaliz- 
able by an orthogonal matrix U{x) with U{x)U'^{x) = 1, 
where 1 denotes the identity matrix. Moreover, the normal- 
ized eigenvectors fi{x) of B{x) form an orthonormal basis, 
fi{x) ■ fj{x) ~ 5ij. Thus, we can construct the trans- 
formation matrix U{x) = {fi{x), . . . , fj^_i{x)) such that 
B{x) = U{x)A{x)U'^{x) where A(a;) is a diagonal matrix 
with the eigenvalues Xi{x) of B{x) along its diagonal. From 
Eq. (HI follows that B{x) is positive definite and hence all 
eigenvalues Xi{x) are positive. Here, we tacitly assume that 
all Tjk are nonzero; if certain transitions are excluded, B{x) 
is positive semidefinite and eigenvalues can be zero. 

Finally, our matrix C{x) is given by 



Cix) =U{x) 
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■U^{x). (6) 
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This standard procedure to diagonalize matrices can be eas- 
ily implemented numerically. However, the diffusion matrix 
B(x) obviously depends on the transition probabilities and 
thus on the abundances of all strategies. Consequently, the 
procedure to calculate C{x) must be continuously repeated as 
time progresses and the state x changes. This is computa- 
tionally inconvenient and therefore, it is desirable to calculate 
C{x) analytically. 

The simplest case with d = 2 strategic types results in a 
one-dimensional Fokker-Planck equation 180. Moreover, in 
special cases, for example in cyclic games such as the sym- 
metric Rock-Scissors-Paper game, the Fokker-Planck equa- 
tion ^ can be approximated in polar coordinates [il3;]. How- 
ever, such cases are non-generic and here we focus on general, 
higher dimensional situations with d > 3. As a particular ex- 
ample to illustrate the framework numerically, we provide a 
detailed analysis of a generic Rock-Scissors-Paper game. 



rf-i 



ik = Ak{x) + \Ckj{x)^j{t) 



(5) 



IV. MULTIPLICATIVE NOISE FOR SPECIFIC 
PROCESSES 



where the ^j (t) represent uncorrected Gaussian white noise 
with unit variance, {^k{t)£,j{t')) = SkjS{t — t'). The matrix 
C{x) is defined by C'^{x)C{x) = B{x) and its off-diagonal 
elements are responsible for correlations in the noise of dif- 
ferent strategic types. 

Fluctuations arising in finite populations are approximated 
by the stochastic term in the Langevin equation Q. For given 



In general, the diffusion matrix B{x) depends not only on 
the frequencies x but also on the payoffs (fitness) of the dif- 
ferent strategic types. In this case, an analytic representation 
of C{x) is only of limited use because it would be valid just 
for one particular game. However, B{x) becomes payoff in- 
dependent if Tj^ (a;) + Tkj{x) is payoff independent. Fortu- 
nately, this holds for the broad and relevant class of pairwise 



comparison processes. In these processes, a focal individual 
/ and a model m with payoffs ttj and tt^ are picked at ran- 
dom and a payoff comparison determines whether the focal 
individual switches its strategy. 

Let j{TTf, 7r„i) be the probability that the focal individual 
adopts the strategy of the model (see e.g. lUJl 1 1511 ) and as- 
sume that every mutation leads to a different strategy. Then 
the transition probabilities from type k to type j read 



^l = 0.0 



^ = 0.2 



Tk]ix) = (1 - ^j.)xkXj-f{T:j,'Kk) + fixk- 



1 



1 



(7) 



for j 7^ k. Consequently, any pairwise comparison process 
with 



7(7rj, TTfe) + 7(7rfc, TTj) = const. 



(8) 



leads to a payoff independent diffusion matrix B{x). If 
Eq. ^ is fulfilled, the noise term in Eq. Q is independent 
of the evolutionary game. In particular, it allows the con- 
sideration of multi-player games in which the payoff func- 
tions are non-linear lll6ll . Examples for evolutionary processes 
that fulfill Eq. (O include the local update process ^ with 
j{ttj , TTfe ) = i +w{ttj ~ TTfc) whcrc w indicates the strength of 
selection acting on payoff differences between different strate- 
gic types. For u; <C 1 selection is weak, payoff differences 
amount to little changes in fitness and the process is domi- 
nated by random updating. For larger w selection strength 
increases but an upper limit is imposed on w by the require- 
ment < 7(7rj, TTfc) < 1. Another example is the Fermi pro- 
cess lfl7| - [l9ll with 7(7r_,, TTfc) = (1 + cxp[-w(7rj - 7rfc)])"\ 
Again w indicates the selection strength but without an upper 
bound. In the limit w — )■ oo, i.e. 7(7rj,7rfc) = [ttj — Tr^] 
where 9 [a;] denotes the Heavyside step function, the imita- 
tion dynamics luOiulll is recovered. A further example, where 
7(7rj, TTfc) does not simply depend on the difference between 
its arguments, is ^{-Kj.-Kk) = 7q^ iHllSl- Incidentally, all 
examples above satisfy 7(7rj,7rfc) + 7(7rfe,7rj) = 1 but, for 
example, there might be resilience to change in the local up- 
date process such that 7(7rj, TTfc) = 0(5+ vj(j^j — TTk)) with 
< a < 1 and hence 7(7rj, ttj.) + 7(7rfe, ttj) = a < 1. 

Last but not least, an example of an important process that 
does not lead to a payoff independent B{x), is given by the 
standard frequency dependent Moran process (J7l |23l |24l or its 
linearized equivalent 1 12, 25]. 

In the following, we concentrate on cases with 7(7rj, TTfc) + 
j{TTk, TTj) = 1. This leads to 



BM^)=j^ 



and 



-XjXk{l - /i) - 



-{Xj +.Tfc) 



for j y^ k 
(9) 



B 



ni^yj^ 



x,{l- x,)(l- A^) + ^(1 + x,{d - 2)) 



Unfortunately, even in this case, a full derivation of C{x) is 
difficult for general d. Thus, we focus on the more manage- 
able but highly illustrative case of d = 3. Henceforth, we set 
xi = X and X2 — y (and x^ — I ~ x — y) for convenience. 
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FIG. 1. (Color online) Elements of the noise matrix C{x, y, z) for 
processes fulfilling Eq. ([8} with d = 3 strategies f or /i = (left) and 
/i > (right), a The element Cxx determines how the noise in the 
x-direction affects the x-coordinate. In the case of /^ = 0, this noise 
vanishes for x — )> 0. For y — i- and 2 — >■ we recover the usual mul- 
tiplicative noise from one-dimensional evolutionary processes Q]. b 
For nonvanishing mutations Cxx increases at x — i- and x- — >■ 1 and 
exceeds the value \— (contour line) in the interior of the simplex, 
which is the maximum value without mutations. The element Cyy 
follows from the transformation x -f^- y. c The element Cxy deter- 
mines how the noise in the x-direction affects the y-coordinate (or 
vice versa). This term is always negative, as explained in the text. 
For X — > or y — !> 0, Cxy vanishes in the case of ^ = 0. For z — !> 0, 



we find Cx 



'Cxx, which ensures that the sum of the noise in 



the x-coordinate, Cxx + Cxy, vanishes on this edge of the simplex: 
If the noise increases x for z = it has to decrease y by the same 
amount, d In the case of /i > 0, the noise no longer vanishes at the 
boundaries of the simplex but the expressions are too cumbersome to 
display explicitly. 



A. No mutations, /i, = 



For d = i and in the absence of mutations, we can give a 
relatively compact analytic expression for C{x), 



B{x) 



1 ( x{l — x) —xy 

N V -^y vi^-y) 



(11) 



(10) xhe eigenvalues of this matrix are 



\±{x) 



K+±L 

2N ' 



(12) 



where K± = x{l — x) ± j/(l — y) and L = \/ K^ + Ax^y'^. 
The eigenvectors are given by 



f±{x) = ^{-K^±L,2xy), 



(13) 



where N = 'JAx^y'^ + {—K- ± L) is a normahzation fac- 
tor From this, it is straightforward to construct the trans- 
formation matrix Z^ (a;) = {f^{x),f_{x)). The product 
U{x)^A{x) l/^ [x) then yields our matrix C{x), which can 
be written as 



^xx KS^ ) 



VQ^ 



VQ] 



i+{^ 



(14a) 



xy 







(14c) 



where Q±± = K± ± L. It is obvious that Cxx{x) > and 
Cyy{x) > 0. Since i > 0, we have Q++ > Q-\ and there- 
fore Cxy{x) < 0. It is remarkable that even for this sim- 
ple B{x), the matrix C{x) already takes a rather complicated 
form that is not easy to interpret. Therefore, the elements 
Cxx{x) and Cxy{x) are plotted in Fig.[T] The remaining el- 
ement, Cyy{x), follows from the symmetry of the matrix, i.e. 
in Fig.[T]a, b the simplex needs to be mirrored along the ver- 
tical axis x ^ y. Recall that the matrix C{x) is independent 
of the evolutionary game that is played and does not depend 
on the microscopic evolutionary process, as long as Eq. (|8), 
j{7Tj,TTk) + j{Trk,TTj) = const., holds. While analytical cal- 
culations of C are also, in principle, feasible for d = 4 and 
d = 5, they lead to very lengthy expressions. 



B. With mutations, /x > 

The procedure to derive C{x) remains the same when in- 
cluding mutations just starting with 



Bix) 



1-A^ fxjl 
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2N 
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-xy 

l+x - 
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-X — y 

i+y 



(15) 



Unfortunately, the analytic expressions for the elements of 
C{x) grow to unwieldy proportions. A graphical illustration 
of the elements of C(a;) for /^ > is shown in Fig. [T] Com- 
pared to (U = 0, the entries of Cxx{x) become larger and the 
noise terms no longer vanish at the boundaries. For example. 



at the corners of the simplex we have for Cxx{x) 



C,,(l,0,0)-^-/ii 



C..(0,1,0) 
C..(0,0,1) 



/To V N 



1 JL 

5V N 



1 J± 



For Cxy{x), we obtain at these points 

Cxy{lAO)^Cxy{0,l,0)^-^^ 
C,y(0,0,l)=0. 

In addition, the functional form of the noise term is altered 
with increasing /i and becomes significantly more complex 
than in the case of no mutations, /i = 0. 

For /i > it is important to mention that for small mutation 
rates serious mathematical intricacies arise in the vicinity of 
absorbing boundaries and saddle-node fixed points ll26ll27ll . 
Intuitively, the reason for these complications arises from the 
diffusion approximation that underlies the derivation of the 
Fokker-Planck equation Q. In systems with such absorb- 
ing boundaries (or sub-spaces where one or more strategic 
types are absent) finite populations spend non-negligible time 
on these boundaries in the limit fj. ^ 1/iV. In contrast, the 
diffusion approximation is based on a continuum approxima- 
tion, which reflects the limit A^ ^ oo, and adds finite size 
corrections for large, finite N. For example, the state of the 
population a; is a continuous variable and hence can be lo- 
cated arbitrarily close to a boundary in both the Fokker-Planck 
or Langevin formalisms, Eq. (|2) and Eq. (|5), but this is im- 
possible in finite populations. More specifically, this implies 
that the case of /i <C 1/N, for which the population could 
get trapped on an absorbing boundary for an extended period 
of time, cannot be captured by the diffusion approximation. 
Consequently, the quality of the approximation is expected to 
decrease for small /i > and to get worse if N is small too. 
Interestingly, this failure of the continuum approximation has 
not been raised in population genetics, where similar consid- 
erations have been made, but typically the description is only 
made on the level of the Fokker-Planck equation and not based 
on stochastic differential equations ||28I1 . However, population 
geneticists are typically interested in simpler scenarios, where 
each type has a fixed fitness. In this case, the corresponding 
deterministic system has no generic trajectories that are close 
to absorbing boundaries. Moreover, the usual diffusion ap- 
proximation in population genetics, where the selection inten- 
sity scales to zero while the population size diverges, leads to 
substantial noise, such that these effects have a minor impact. 

For evolutionary games with small mutation rates one is 
typically not forced to apply this relatively complex approxi- 
mation. Instead, a time scale separation between mutation and 
selection allows to approximate the dynamics based on a pair- 
wise consideration of strategies by considering an embedded 
Markov chain over the (quasi absorbing) homogenous states 
of the population ll29l - [33ll . 
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V. APPLICATION TO THE ROCK-PAPER-SCISSORS 
GAME 

To compare our approach based on stochastic differential 
equations in d — 1 dimensions to individual based simulations 
with d strategies, we focus on the case of d = 3 where we 
have obtained closed analytical results above for /i = 0. As 
an example, we consider the cyclic dynamics in the Rock- 
Scissors-Paper-Game, which is not o nly a po pular children's 
game, but also relevant in biological Il34l437ll and social sys- 
tems Il38ll39ll . Moreover, the evolutionary dynamics of this 
game is theoretically very well understood, both in infinite as 
well as in finite populations Il2l [l3l[l4ll40l - l42ll . Here we focus 
on a generic Rock-Scissors-Paper game with payoff matrix 



(16) 



which avoids artificial symmetries. According to the repli- 
cator equation, the game exhibits saddle node fixed points at 
X = l,y = 1, and ^ = 1 — x — y = las well as an interior 
fixed point at a; = {\,\,\), independent of the parameter s 
but s controls the stability of x. For s > 1, a; is a stable fo- 
cus and an unstable focus for s < 1. We do not consider the 
non-generic case of s = 1, which exhibits closed orbits pt). 

Let us first analyze the fixation probabilities and the fixa- 
tion times for the case without mutations /^ = for a game 
with s ~ 1.4, such that x is an attractor The replicator equa- 
tion, obtained in the limit iV ^- oo, predicts that fixation never 
occurs and that the population evolves towards x, see Fig.|2^. 
However, in the stochastic system ultimately two strategies are 
lost, see Fig.llJ), c. Most importantly, the sample trajectories 
generated by the stochastic differential equation (|5]l in Fig.|2]3 
do not exhibit any qualitative differences when compared to 
individual based simulations. Fig. |2];. As a complementary 
scenario, we consider a game with s = 0.8, such that a; is a 
repellor According to the replicator equation trajectories now 
spiral away from x towards the boundaries of the simplex and 
approach a heteroclinic cycle. The probabilities that the pop- 
ulation reaches any one of the three homogenous absorbing 
states in either scenario is shown in Fig. |3] together with the 
average time to fixation. 

In all cases excellent agreement between the Langevin 
framework (|5j and individual based simulations is obtained. 
Note that larger A^ not only improve the approximation Eq. (|5) 
but also result in a performance gain as compared to individ- 
ual based simulations. More specifically, the computational 
effort scales linearly with population size for simulations but 
remains constant when integrating Eq. (|5]l numerically. 

For non-vanishing mutation rates, /i > 0, absorption is no 
longer possible and the boundary of the simplex becomes re- 
pelling. Fig. m depicts the average frequencies of each strate- 
gic type as a function of yu. For small /i the population can 
still get trapped for considerable time along the boundary, 
which results in systematic deviations between the stochas- 
tic differential equation (|5]l and individual based simulations. 




FIG. 2. Comparison of trajectories for the a deterministic case 
(replicator equation), b stochastic differential equation, Eq. ([5}, and 
c individual based simulations for an RSP interaction with weakly 
attracting interior fixed point x with s — 1.4 (c.f. Eq. ( |16b ) and 
A^ — 1000 (in b, c). Stochastic trajectories start close to x and end 
at time T = 349.9 (b) and T = 395.8 (c). Simulations can also be 
performed online at 114 311 . 



Fig. |4] More specifically, for n < 1/N the deviations in- 
crease with decreasing fj, but the agreement remains excellent 
for /i > 1/A^. This threshold simply means that for larger /i 
the time spent along the boundary can be neglected. In the 
limit N ^ oo the replicator-mutator Eq. (|5) exhibits a stable 
limit cycle H . 



VI. DISCUSSION 

Evolutionary dynamics can be implemented in multiple 
ways. In particular, if only two types are present, the dynamics 
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FIG. 3. Probability of absorption and time to absorption (in the 
absence of mutations, /i = 0). a, c Probability to reach one of the 
three absorbing homogenous states of rock (solid line, ■), scissors 
(dashed line. A), or paper (dotted line, ♦) as well as b, d the as- 
sociated time to absorption for rock as a function of the population 
size A^ with a stable (top row, s — 1.4) and unstable (bottom row, 
s = 0.8) interior fixed point x. If x is an attractor (a, b) it may take 
exceedingly long times for large A^ until the population reaches an 
absorbing state - even though this will inevitably occur. Because of 
this results are only shown up to A'' = 3000. No such limitations 
occur if a; is a repellor. The symbols indicate results from individ- 
ual based simulations. Error bars and grey shaded areas indicate the 
standard deviation of the mean for simulations and stochastic cal- 
culations, respectively. Simulation results and results based on the 
stochastic differential equation Eq. ^ were both averaged over 10'' 
independent runs. 



reduces to a single dimension and is typically solvable analyt- 
ically, even if the population is finite and demographic noise 
is present. If more than two types are present, it is signifi- 
cantly more challenging to describe the dynamics analytically. 
When the mutation rates are sufficiently small OSll . there are 
typically at most two types present at any time, thus leading to 
situations which justify simpler tools based on the interaction 
of two types. In the limit of infinite populations, determinis- 
tic differential equations arise and can be used to describe the 
system even when the population size is large, but finite. How- 
ever, it is more challenging to incorporate noise in this case. 
Here, we have proposed a way to address this issue by deriv- 
ing a Langevin equation for more than two types. For large 
populations, it is typically much more efficient to solve these 
equations numerically than to resort to individual based sim- 
ulations. A detailed performance comparison is shown in Fig. 
|5] Computational costs of simulations increase slowly with 
the number of strategic types, ~ d^/^, but increase with the 
population size as N'^. In contrast, SDK's are essentially un- 
affected by N, but computational costs arise from repeatedly 
solving for eigenvalues and eigenvectors of B{x). In theory, 
analytical solutions are available for d < 5 but are probably 
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FIG. 4. Average frequency (solid line) and standard deviation (grey 
shaded area) of the three strategies rock (a, b), scissors (c, d), and 
paper (e, f) as a function of the mutation rate p for small (left col- 
umn, A^ — 100) and large populations (right column. A' = 10^) 
with a stable interior fixed point x, s — 1.4, based on Eq. (|5]l. For 
comparison, symbols and eiTor bars depict results from individual 
based simulations. No simulation results are shown for A'^ = 10'' be- 
cause of prohibitive computational efforts. Simulations and Eq. (|5} 
are in excellent agreement for p > 1/N but for smaller p substan- 
tial deviations occur (see text for details). Stochastic fluctuations 
decrease with increasing N and for A'^ = lO"" the population spends 
most of the time in the close vicinity of x. Simulation results and re- 
sults based on the stochastic differential equation Eq. (|5) were both 
averaged over 10^" time steps after a relaxation time of lO'' steps 
(dt = 0.01 for the Langevin equation). 



meaningful in practice only for d = 2, 3 and numerical meth- 
ods are required for d > 3, which scale with d'^ HS\ . To 
illustrate this, for the data point N = 10000 in Fig. [3]:, d the 
simulations required approximately two months (1418 hours) 
to complete as compared to 49 minutes for the corresponding 
calculation based on stochastic differential equations, which 
corresponds to a 1700-fold performance gain. However, if 
only a small number of one type is present in a population our 
approach does not work well unless there are no mutations 
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or the product of the population size and the mutation rate is 
sufficiently high, N n > 1, such that the time spent along the 
boundary becomes negligible, cf. Fig.|4] 

In summary, our approach establishes a transparent link be- 
tween deterministic models and individual-based simulations 
of evolutionary processes. Moreover, the resulting stochastic 
differential equations provide substantial speed-up compared 
to simulations, and therefore may serve well in the investiga- 
tion of multidimensional evolutionary dynamics. 
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FIG. 5. (color online) Performance comparison of individual based 
simulations (IBS) versus stochastic differential equations (SDE). a 
ratio of the CPU times CPUsde/CPUibs as a function of the popu- 
lation size, A'^, and the number of strategic types, d. The bold contour 
indicates equal performance. For small A'^ and large d IBS are faster 
(left of solid line), but for larger TV and smaller d SDE are faster 
(right of solid line). Each contour indicates a performance difference 
of one order of magnitude. Note that IBS for d = 2, 3 are based 
on analytical calculations of the eigenvalues/eigenvectors of B{x), 
which requires substantially less time and hence explains the dis- 
continuity in the contours, b computational time with d = 10 as a 
function of N for IBS and SDE. As a reference for the scaling A'^^ 
and a constant are shown, c computational time with A'^ = 5000 
as a function of d for IBS and SDE. As a reference for the scaling 
d^'^ and d"^ are shown. For a proper scaling argument much larger 
d are required but already d = 100 far exceeds typical evolutionary 
models and hence is only of limited relevance in the cuiTent context. 
All comparisons use a constant payoff matrix and the local update 
process (SD (such that Ak (x) = and 7(7rj ,TVk)~ 1/2), a mutation 
rate of 1/A' and are based on at least 1000 time steps as well as at 
least one minute running time. CPU time is measured in millisec- 
onds required to calculate 1000 time steps. The time increment for 
the SDE is dt = 0.01. 
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